g = 0;  % Polarizer角度
f = 90; % Retarder相位延迟

t1_range = linspace(0, 360, 100); % t1从0到360度
t2_range = linspace(0, 360, 100); % t2从0到360度

S0_output = zeros(length(t1_range), length(t2_range));

for i = 1:length(t1_range)
    for j = 1:length(t2_range)
        t1 = t1_range(i);
        t2 = t2_range(j);

        rotated_polarizer = 0.5 * [1, cosd(2*g)*cosd(2*t1), cosd(2*g)*sind(2*t1), 0;
            cosd(2*g)*cosd(2*t1), cosd(2*t1)^2 + sind(2*g)*sind(2*t1)^2, (1 - sind(2*g))*sind(2*t1)*cosd(2*t1), 0;
            cosd(2*g)*sind(2*t1), (1 - sind(2*g))*sind(2*t1)*cosd(2*t1), sind(2*t1)^2 + sind(2*g)*cosd(2*t1)^2, 0;
            0, 0, 0, sind(2*g)];
        
        rotated_retarder = [1, 0, 0, 0;
            0, (cosd(2*t2))^2 + cosd(f)*sind(2*t2)^2, (1 - cosd(f))*sind(2*t2)*cosd(2*t2), -sind(f)*sind(2*t2);
            0, (1 - cosd(f))*sind(2*t2)*cosd(2*t2), sind(2*t2)^2 + cosd(f)*cosd(2*t2)^2, sind(f)*cosd(2*t2);
            0, sind(f)*sind(2*t2), -sind(f)*cosd(2*t2), cosd(f)];

        S0 = [1; 1; 0; 0];

        S_out = rotated_polarizer * (rotated_retarder * S0);

        S0_output(i, j) = S_out(1);
    end
end

figure;
mesh(t1_range, t2_range, S0_output');
xlabel('t1 (degrees)');
ylabel('t2 (degrees)');
zlabel('Emerging S0 Intensity');
title('S0 Intensity Variation for 1/4 retarder');
